Nonlinear phase mixing and phase-space cascade of entropy in gyrokinetic plasma turbulence 
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Electrostatic turbulence in weakly collisional, magnetized plasma can be interpreted as a cascade of entropy 
in phase space, which is proposed as a universal mechanism for dissipation of energy in magnetized plasma 
turbulence. When the nonlinear decorrelation time at the scale of the thermal Larmor radius is shorter than 
the collision time, a broad spectrum of fluctuations at sub-Larmor scales is numerically found in velocity and 
position space, with theoretically predicted scalings. The results are important because they identify what is 
probably a universal Kolmogorov-like regime for kinetic turbulence; and because any physical process that 
produces fluctuations of the gyrophase-independent part of the distribution function may, via the entropy cas- 
cade, result in turbulent heating at a rate that increases with the fluctuation amplitude, but is independent of the 
collision frequency. 

PACS numbers: 52.30.Gz, 52.35.Ra, 52.65.Tt 



43 

Oh- 



4h 1 
O 1 

43 : 



> 

OO 

in 

(N 



oo 
o 



X 



Introduction. Turbulence is inherently nonlinear and dy- 
namically complicated. In the general case, a broad spectrum 
of fluctuations is excited, in both wave number and frequency. 
For turbulent, magnetized plasma, the equations of magne- 
tohydrodynamics provide a pedagogically rich description of 
the dynamics. However, for those turbulent eddies whose par- 
allel wavelengths (relative to the magnetic field) are compara- 
ble to or smaller than the collisional mean free path and whose 
perpendicular wavelengths are comparable to or smaller than 
the Larmor radius of one of the constituent species of the 
plasma, magnetohydrodynamic theory breaks down. In such 
cases, the gyrokinetic (GK) theory (H 0] represents a rigor- 
ous limit of plasma kinetics for anisotropic (Jk\\ <C k±), low- 
frequency (u) <C f2, the ion cyclotron frequency) fluctuations. 
In this Letter, we present a GK description of turbulence in 
a simplified situation, chosen to isolate a novel phenomenon 
which is a generic component of all GK turbulence: the si- 
multaneous cascade of entropy to smaller scales in both real 
space and velocity space. This phase-space cascade is the 
mechanism by which turbulent energy associated with fluctu- 
ating fields is brought to small scales in velocity space, where 
even very infrequent collisions are sufficient to provide irre- 
versibility and thus heating. Below, we present the theory and 
first-principles simulations of the phase-space cascade in a ho- 
mogeneous, electrostatic, magnetized plasma. 

It is well known that Landau and Barnes damping of elec- 
tromagnetic plasma fluctuations lead to the generation of 
small-scale structures in f(v»), where / is the one-particle 
distribution function, and uii is the velocity coordinate along 
the background magnetic field Q 0). This is associated 
with the free-streaming of particles along the field. As t in- 
creases, a single Fourier harmonic of the distribution func- 
tion fk n ~ e n^n* gets progressively more oscillatory in v\\- 
space. Eventually, even infrequent collisions are sufficient to 
smooth these oscillatory features, since the collision opera- 
tor is roughly a diffusion operator in velocity space. As long 
as collisions are sufficiently infrequent, the damping rate de- 



pends not on the collision rate, but on the nature of the wave 
and its phase velocity relative to the thermal speeds of the 
plasma species. Physically, Landau damping is the smearing 
of spatial perturbations that occurs when there is a spread in 
the distribution of parallel velocities. We recall for future ref- 
erence that this generation of velocity-space structure is inde- 
pendent of the fluctuation amplitudes. 

Besides this linear parallel phase mixing, there exists a non- 
linear phase mixing process ^ that, in a strongly turbulent 
plasma and at spatial scales smaller than the Larmor radius, 
drives the formation of structure in f(v±) much more rapidly 
than parallel phase mixing drives f(v\\). Physically, this non- 
linear phase mixing is the smearing of spatial perturbations 
due to the spread in the distribution of gyroaveraged E x B 
velocities (see Fig. []]). Unlike for the parallel phase mixing, 
the rate of generation of v-space structure by this process is 
proportional to the fluctuation amplitude. In this Letter, we 
present a study of this nonlinear process, which we interpret as 
a turbulent cascade of entropy in phase space |6i]. As such, it 
represents a conceptually novel nonlinear phenomenon, where 
generation of small scales in the position and velocity space 
occurs in an intertwined way. This process, which is likely 
to be a fundamental and ubiquitous feature of magnetized 
plasma turbulence, has never been numerically diagnosed and 
analyzed before, although Krommes |3] did point out the gen- 
eral possibility of the coupling between position and velocity 
space. 

Gyrokinetics in 2D. Let the distribution function be / = 
Fq + Sf, where Fo is a Maxwellian with density no and tem- 
perature To, and Sf = h — qtpFo/To, where q is the particle 
charge and <p is the electrostatic potential. To keep the focus 
on the nonlinear process, we consider electrostatic GK turbu- 
lence in slab geometry with ku = 0. Then the non-Boltzmann 
part h of the perturbed ion distribution function satisfies (H 
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FIG. 1 : Schematic view of the nonlinear phase mixing superimposed 
on the potential from the Run (iii) (see Table[l]l at i/rinit = 10 and 
the largest wavelength mode taken out. When the fluctuation scale 
I < p, the gyroaverage of the electric field induces a decorrelation 
of the distribution function at the velocity-space scale corresponding 
to the difference in Larmor radii £ v = 8v/£l ~ £ [see @], 



where Bq is the background magnetic field aligned with the 
z-axis and (•) r is the gyroaverage holding the guiding center 
position R constant. The collision operator C[h] used in our 
simulations contains pitch-angle scatteringand energy diffu- 
sion with proper conservation properties |7||. The quasineu- 
trality condition yields 



q J (h) r dv = qY,e lk r J Jo 



fi 



hk dv, 



(2) 

where (-) r denotes the gyroaverage at fixed particle posi- 
tion r, Jo is the Bessel function, Q = ^ s q 2 s no s /T 0s for 
Boltzmann-response (3D) electrons or Q = q 2 noi/Toi for 
no-response (2D) electrons, and s and i are the species in- 
dices. Our results are not affected by the choice of the electron 
response. For concreteness, we henceforth use no-response 
electrons since electrons cannot contribute to the potential if 
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exactly. In the absence of collisions, the system has 



two positive definite conserved integrals 
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where T = Io{k\p 2 j2)e~ k ^ p ' 1 1 2 , I is the modified Bessel 
function and p is the ion thermal Larmor radius. The invariant 
38 is proportional to minus theperturbed part of the entropy of 
the system, — j / In / dr dv yl |4[] . Here we will refer to 38 as 
"entropy" to emphasize this connection. The second invariant 
€ is conserved in the 2D electrostatic case only. 

Scalings. A scaling theory of the entropy cascade in the 
sub-Larmor scale range can be developed in a way reminiscent 
of the Kolmogorov-style turbulence theories J(|. Assume that 



at (perpendicular) scales £ <C p, the transfer of entropy is local 
in scale. On dimensional grounds, the entropy flux is 
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until it reaches the collisional dissipation scale, where vth is 
the thermal speed and Tg is the nonlinear decorrelation time at 
scale I. The neglect of the ip 2 term in 38 [see (O] is justified 
post hoc due to its smallness in the I <C p regime [see (O and 
Fig- Eta)]. There is a self-consistent electrostatic potential at 
the scale I: from (O, 
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Here we have assumed that the nonlinear phase mixing pro- 
duces velocity-space structures correlated with the spatial 
scale via (see Fig.[T]and Refs. |0,[3l) 



Vth 



(7) 



This has allowed us to estimate the velocity integral in © as a 
random-walk-like accumulation of the integrand represented 
by the product of hk, which is random function of v± whose 
"step size" is given by (0 with I ~ fcj , and of the Bessel 
function, which introduces a reduction factor of (i/p) 1 ' 2 . 

The decorrelation time ti may be estimated by balancing 
the dt term with the nonlinear term in (Q]i, leading to 
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Substituting © into © and © into © yields h 
ip ~ £ 7 / 6 . Therefore, the spectra of h and <p are 
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where 



k± jT \h k \ 2 /2F dv and 



v ( k ±) = J2\ k± \=k ± q 2 no\tp k \ 2 /2To. Note that the total en- 
tropy (O can be expressed as 38 = f[Ef l (k±) — E v (kj_)] dk±. 

Dissipation Cutoff. From the balance between the nonlin- 
ear decorrelation time ((HJl and the collision time v^ 1 , one ob- 
tains an estimate of the dissipation cutoff scales in both the 
velocity space and the real space [see (0)]. Using C[h] ~ 
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where t p is the nonlinear decorrelation time measured at 
I = p. We have introduced a new dimensionless number 
D to characterize the scale separation in gyrokinetic turbu- 
lence: analogous to the Reynolds number in fluid turbulence, 
large D corresponds to a broader scaling range over which the 
entropy cascade extends, and to dissipation at smaller scales. 
Here, however, the smallest spatial scale observed is deter- 
mined by the u-space scale for which diffusion in velocities 
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TABLE I: Index of the runs. 



Run N x x N y 


N, x 2iV A 




D 


k± c p 


(i) 64 2 


32 2 


5.6- 10" 3 


48 


20 


(ii) 128 2 


64 2 


1.9 ■ 1(T 3 


118 


35 


(iii) 256 2 


128 2 


7.4 ■ 1(T 4 


440 


77 
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FIG. 2: Time evolution of 38 and € [Eqs. l[3} and lO] normalized to 
initial 38. The runs (i)-(iii) are indexed in Table U Evolution of € 
does not differ among runs significantly, and is given for run (iii). 

becomes important, through the correlation between real and 
velocity space given by (01. The fact that D increases with the 
amplitude of the fluctuations at the Larmor scale clearly dis- 
tinguishes this process from linear Landau damping. We note 
that for 3D gyrokinetic turbulence, the nonlinear phase mixing 
is a much faster process than the linear one if the fluctuation 
amplitude is sufficiently large. 

Numerical Simulations. We now report the first-of-a-kind 
numerical investigation of the entropy cascade in phase space, 
carried out with the GK code AstroGK. The code uses a 
Fourier pseudo-spectral scheme for the real-space dimensions 
perpendicular to the background magnetic field and a Leg- 
endre collocation scheme for the velocity-space integrations. 
The velocity space is discretized in energy e = v 2 and 
A = v\j e. In the absence of collisions, AstroGK conserves 
the invariants (0 and (0) with a high precision. 

The results reported below were obtained in three runs at 
decreasing collision frequency v and correspondingly increas- 
ing spatial and velocity resolution. They are indexed in Ta- 
ble U where N x x N y is number of collocation points in 
the real space and N £ x 2N\ is the number of grid points 
in velocity space — the factor of 2 corresponds to the sign 
of uii = ±y/s(l — A). Our highest-resolved run required 36 
wallclock hours on 8192 processors. 

The code evolves g = h — qFo((p)n/To and ip via Eqs. (Q]i 
and 10. We take the box size L x = L y = 2irp and start 
from the initial condition ginit = ga[cos(2x / p) + cos(2y / p) + 
x(x,u)]Fq, where go is a constant and x(a;,y) is a small- 
amplitude white noise superimposed on all Fourier modes. 



From (0, we can calculate tpinit- 

Time Evolution. The initial \k x p\, \k y p\ = 2 configu- 
ration is unstable: the amplitudes of <p corresponding to 
\k x p\, \k y p\ = 1 grow and then saturate around t/T- m n, ~ 
9, where r init = 2kB / (ck\\\{wivit)ii\\) is the turnover 
time associated with the initial condition and ||(<^)h|| = 
[(1/nrj) // \{<*p)r\ 2 Fq dv dR] 1 / 2 . The nonlinear interactions 
between modes produce smaller scales down to a cutoff de- 
termined by D [see (TTObl . The turbulent spectra fill up by 
t/r in it ~ 10, then decay with time. 

The time evolution of the collisionless conserved quantities 
38 and € [see (01 and (0)] is shown in Fig. [2] During the initial 
growth of the instability of the \k x p\, \k y p\ = 1, 38 decays 
very slowly at a rate ~ u, consistent with a collisional decay 
rate associated with the large-scale phase-space variation of 
^init- Once turbulence develops, 38 decays more rapidly as 
the entropy cascade transfers it nonlinearly to smaller scales in 
phase space, until the fluctuations of the distribution function 
are thermalized (dissipated) at the collisional cutoff. 

The decrease of 38 from its initial value corresponds to the 
amount of entropy (heat) production due to the irreversible 
collisional smearing of the distribution function. The turbu- 
lence that follows the initial instability enhances the heating, 
suggesting that small-scale velocity-space structure is gener- 
ated (this is confirmed below). As expected, the rate of dissi- 
pation is not strongly affected by the collision frequency, i.e., 
there is a finite amount of dissipation even as the collision 
frequency tends to zero. The dissipation rate is determined 
instead by the nonlinear cascade rate. 

While 38 decays, € stays almost constant. If we increase the 
size of the simulation box, the \k x p\, \k v p\ = 1 modes them- 
selves become unstable to even longer-wavelength modes. We 
attribute both this instability and the failure of € to decay to 
the intrinsic tendency of € to have an inverse cascade 
which we do not discuss here. 

Spectra and Scalings. The wave-number spectra of the 
decaying developed turbulence are given in Fig. 0a). They 
are angle integrated over wave-number shells |fej_| = k±, nor- 
malized by 38(i) at each time and then averaged over time for 
10 < t/ri n it < 15. As resolution is increased, the spectra 
appear to converge to the theoretically predicted scalings (O, 
which supports the validity of our dimensional and physical 
considerations of the entropy cascade. 

To characterize the entropy cascade in the velocity space, 
Plunk et al. JH introduced velocity-space spectra E g (p) = 
J2kP\9k(p)\ 2 , where g k (p) = J Ja(pv±)g k (v) dv is a Han- 
kel transform. The theoretical expectation is that E g (p) ~ 
p- 4 /3 because the real- and velocity-space scales should be 
related according to <J7J, which, in terms of the dual variable 
p, becomes k±p ~ pvth- The time-averaged Hankel spectrum 
E g {p) obtained in our simulations is shown in Fig. 0b). This 
again shows approximate consistency with the theoretical pre- 
diction and confirms that small-scale structure is formed in the 
velocity space. 

Dissipation Cutoff. In Table [I] we show for each of our 
runs the dimensionless number D = (vTp)^ 1 , where r p = 
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FIG. 3: Time-averaged normalized (a) wave-number (Fourier) spec- 
tra E h (kx)/$B and E v (k±)/%& [cf. ©], and (b) velocity-space 
(Hankel) spectrum E g (p) /SB for the runs indexed in Table U Theo- 
retically predicted slopes are given for comparison. 



27t.Bo/(c/s'^||(i/2 / }r||) measured at t/r m it = 10 and ip' is 
ip with the \k x p\, \k y p\ = 1 modes taken out (see also Fig. 
Q}. Also shown is the theoretical estimate (TTOb for the wave- 
number cutoff k± c p = otD z f b , where a = 2 is an empirical 
value that corresponds to our particular set up. Comparing 
with the wave-number and velocity-space spectra in Fig. [5J 
we see that ( fTOb describes the resolution requirements quite 
well. With fewer velocity grid points, we find shallower wave- 
number spectra than the resolved ones, while with more, we 
resolve below the velocity cutoff without any change in the 
wave-number spectra. Thus D is a good indicator of neces- 
sary and sufficient resolution in full 4D phase space. 

Conclusions. We have presented electrostatic, decaying 
turbulence simulations for weakly collisional, magnetized 
plasmas using the gyrokinetic model in 4D phase space 
(two real-space and two velocity-space dimensions). Landau 
damping was removed from the system by ignoring variation 
along the background magnetic field. Nonlinear interactions 
introduce an amplitude-dependent perpendicular phase mix- 
ing of the gyrophase-independent part of the perturbed distri- 
bution function and create structure in v i which is finer for 



higher k±. We have found that the wave-number (Fourier) 
and velocity-space (Hankel) spectra of the perturbed distribu- 
tion function and the resulting electrostatic fluctuations at sub- 
Larmor scales agree well with theoretical predictions based on 
the interpretation of the nonlinear phase mixing as a cascade 
of entropy in phase space @] ■ We have introduced a di- 
mensionless number D (analogous to Reynolds number) that 
characterizes the scale separation between the thermal Larmor 
scale and the collisional cutoff in phase space [see dlOH . and 
showed that this number correctly predicts the resolution re- 
quirements for our simulations. 

We note that there are, in general, entropy cascades for each 
plasma species. Equations for the gyrokinetic turbulence at 
and below the electron Larmor scale are mathematically simi- 
lar to the model simulated here and identical arguments apply 
JfJ, @]. Similar considerations are also possible for ion-scale 
electromagnetic turbulence j^l and for minority species. 

The small-scale phase-space structure that we have discov- 
ered is likely to be a universal feature of strong, magnetized 
plasma turbulence. Understanding it theoretically and diag- 
nosing it numerically is akin to the inertial-range studies for 
Kolmogorov turbulence, extended to the kinetic phase space. 
One should expect rich and interesting physics to emerge and 
it is likely that predicting large-scale dynamics will require 
effective models for the small-scale cascade. An immediate 
key physical implication of the existence of the entropy cas- 
cade is a turbulent heating rate independent of collisionality 
in weakly collisional plasmas. 

Numerical computations were performed at NERSC, 
NCSA and TACC. This work was supported by the U.S. DOE 
Center for Multiscale Plasma Dynamics, STFC, and the Lev- 
erhulme Trust Network for Magnetised Plasma Turbulence. 
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